home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Correlation.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-10-18  |  14.4 KB  |  774 lines
  1. /* Correlation Statistics */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30) 
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.  /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols>2 then do
  43.     'Message "Only two columns allowed!"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47.  
  48. /* Get cell reference for output range */
  49. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  50. if out_cell="" then do
  51.     'DEFPUBSCREEN("Workbench")'
  52.     exit
  53. end
  54. if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
  55.     'Message "Invalid cell reference"'
  56.     'DEFPUBSCREEN("Workbench")'
  57.     exit
  58. end
  59. /* Suppress Screen Redraw to Speed Things Up */
  60. 'Refresh 0'
  61.  
  62. /* Open a small output window on tcalc screen*/
  63. fo=0
  64. CR='0a'x
  65. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  66. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  67.      call writeln(6Info, DisplayMsg)
  68.     fo=1
  69. end
  70. else do
  71.     'Message "TCALC Screen not available for Progress messages"'
  72. end
  73. CALL DELAY(150)
  74.  
  75. /* Get cell references for top cell in each column */
  76. 'SelectCell' start_cell
  77. do col=start_col to end_col
  78.     'GetCursorPos'
  79.     top_cell.col=result
  80.     'Column 1'
  81. end
  82. if fo then call writech(6Info,"Progress...0 ")
  83. /* Get labels for later use on output */
  84. 'SelectCell' start_cell
  85. 'GetValue'
  86. testlabel=result
  87. testlabel=strip(testlabel)
  88. if datatype(testlabel,'n')=1 then do
  89.     labelflag=0
  90.     do x=1 to NCols
  91.     title.x="Column "||x
  92.     end
  93. end
  94. else do
  95.     labelflag=1
  96.     NRows=NRows-1
  97.     do x=1 to NCols
  98.     'GetValue'
  99.     str=result
  100.     title.x=translate(strip(str),"_"," ")
  101.     'Column 1'
  102.     end
  103. end
  104. if fo then call writech(6Info,"10 ")
  105.  
  106. /* Get data from cell range */
  107. col=start_col
  108. lav=0
  109. tot=0
  110. count.=0
  111. total.=0
  112. do x=1 to NCols
  113.     'SelectCell' top_cell.col
  114.     if labelflag=1 then 'CursorDown 1'
  115.     do y=1 to NRows
  116.         'GetValue'
  117.         valtest=result
  118.         if datatype(valtest)='NUM' then do
  119.             'GetValue'
  120.             val=result
  121.             data.x.y=val
  122.             tot=tot+val
  123.             total.x=tot
  124.             count.x=1+count.x
  125.         end
  126.         'CursorDown 1'
  127.     end
  128.     col=col+1
  129.     tot=0
  130.     lav=0
  131.     val=0
  132. end
  133. if fo then call writech(6Info,"20 ")
  134. if count.1 ~= count.2 then do
  135.     'Message "Equal sample sizes required - aborting!"'
  136.     'DEFPUBSCREEN("Workbench")'
  137.     exit
  138. end
  139.  
  140. /* Calculate Means */
  141. mean.=0
  142. do x=1 to NCols
  143.     mean.x=total.x/count.x
  144. end
  145. if fo then call writech(6Info,"30 ")
  146.  
  147. /* Calculate Standard deviation and Variance */
  148. dat=0
  149. meenx=0
  150. tot.=0
  151. Vsq.=0
  152. diff.=0 /* Array holding difference between value and column mean */
  153. sum.=0 /* Array holding sum of value minus mean of value squared */
  154. sd.=0 /* Standard deviation array */
  155. var.=0 /* Variance array */
  156. do x=1 to NCols
  157.     sum.x=0
  158.     meenx=mean.x
  159.     do y =1 to Nrows
  160.         dat=data.x.y
  161.         diff.x.y=dat-meenx
  162.         sum.x=(dat-meenx)**2+(sum.x)
  163.         tot.x=tot.x+data.x.y
  164.         Vsq.x=Vsq.x+(data.x.y)**2
  165.     end
  166.     N=(count.x)-1
  167.     var.x=(sum.x)/N
  168.     sd.x=sqrt(var.x)
  169. end
  170. if fo then call writech(6Info,"40 ")
  171. N=NRows
  172. xy=0
  173. summulti=0 /* Sum of product of diffs */
  174. diffMult.=0 /* Array holding values of diffx times diffy */
  175. do y=1 to N
  176.     diffMult.y=(diff.1.y)*(diff.2.y)
  177.     summulti=summulti+(diffMult.y)
  178.     xy=xy+(data.1.y)*(data.2.y)
  179. end
  180.  
  181. if fo then call writech(6Info,"50 ")
  182.  
  183.  
  184. /* Calculate Pearson r */
  185. rs=0
  186. pr=0
  187. pr=summulti/sqrt((sum.1)*(sum.2))
  188. rs=((summulti)**2)/((sum.1)*(sum.2))
  189. pr=sqrt(rs)
  190. /*top=(N*xy)-(tot.1)*(tot.2)
  191. bot=sqrt((N*(Vsq.1)-(tot.1)**2)*(N*(Vsq.2)-(tot.2)**2))
  192. pr=top/bot
  193. rs=pr**2*/
  194.  
  195. /* Calculate Confidence Interval for Pearson R */
  196.  
  197. Zr=.5*LOG(1+pr)-.5*LOG(1-pr)
  198. JJ=sqrt(1/(n-3))
  199. CIP95=Zr+1.96*JJ
  200. CIN95=Zr-1.96*JJ
  201. CIP99=Zr+2.58*JJ
  202. CIN99=Zr-2.58*JJ
  203. if fo then call writech(6Info,"60 ")
  204.  
  205. /* Test Against True Population Correlation */
  206.  
  207. PCRR=rtgetstring("0","Enter True Population Correlation (Optional)","Input Request")
  208. PCORR=strip(PCRR)
  209. if PCORR>.99 then PCORR=.99
  210. if rtresult~=0 & PCORR~=0 then Do
  211.     Zp=.5*LOG(1+PCORR)-.5*LOG(1-PCORR)
  212.     Zed=(Zr-Zp)/JJ
  213.     Call NPROB(zed)
  214. End
  215.  
  216. /* Calculate standard Error */
  217. sterr=0
  218. t=0
  219. df=N-2 /* Degrees of Freedom */
  220. sterr=sqrt((1-rs)/df)
  221. t=pr/sterr
  222.  
  223. /* Calculate T distribution function */
  224.  
  225. P=PROBT(t,df)
  226. if t>=0 then P=1-P
  227. PT=P*2
  228. /* Calculate T Critical */
  229.  
  230. TCRIT1=TCRITICAL(.95,df)
  231. TCRIT2=TCRITICAL(.99,df)
  232. TCRIT3=TCRITICAL(.975,df)
  233. if fo then call writech(6Info,"70 ")
  234. TCRIT4=TCRITICAL(.995,df)
  235.  
  236. /* Calculate Column Squares and Products */
  237. SQ.=0
  238. PR.=0
  239. TSQ.=0
  240. TPR.=0
  241. do x=1 to NCols
  242.     do y=1 to NRows
  243.         SQ.x=(data.x.y)**2
  244.         TSQ.x=(TSQ.x)+(SQ.x)
  245.     end
  246. end
  247. if fo then call writech(6Info,"80 ")
  248. do x=1 to NCols
  249.     do z=1 to NCols
  250.         do y=1 to NRows
  251.             PR.x.z=(data.x.y)*(data.z.y)
  252.             TPR.x.z=(TPR.x.z)+(PR.x.z)
  253.         end
  254.     end
  255. end
  256. N=count.1
  257. if fo then call writech(6Info,"90 ")
  258. /* Calculate Correlation */
  259. top.=0
  260. bot.=0
  261. corr.=0
  262. do x =1 to NCols
  263.     do z=1 to NCols
  264.     top.x.z=(N*(TPR.x.z))-((total.x)*(total.z))
  265.     bot.x.z=sqrt((N*(TSQ.x)-(total.x)**2)*(N*(TSQ.z)-(total.z)**2))
  266.     corr.x.z=(top.x.z)/(bot.x.z)
  267.     end
  268. end
  269. if fo then do
  270.     call writeln(6Info,"100 ")
  271.     call writeln(6Info,"Writing output to window...")
  272. end
  273.  
  274. /* Output */
  275. 'SelectCell' out_cell
  276. 'ColumnWidth' 15
  277. 'Put "Correlation: UnGrouped Data"'
  278. 'CursorDown 2'
  279. 'Column 1'
  280. do x=1 to NCols
  281.     'GetCursorPos'
  282.     first_cell.x=result
  283.     title=""""||title.x||""""
  284.     'Alignment 2'
  285.     'Put' title
  286.     'Column 1'
  287. end
  288. 'SelectCell' first_cell.1
  289. 'Column -1'
  290. 'CursorDown 1'
  291. do x=1 to NCols
  292.     title=""""||title.x||""""
  293.     'Alignment 2'
  294.     'Put' title
  295.     'CursorDown 1'
  296. end
  297. d=0
  298. b=0
  299. do x=1 to NCols
  300. 'SelectCell' first_cell.x
  301. 'CursorDown 1'
  302.     do z=1 to NCols
  303.         'Put' format(corr.x.z,,4)
  304.         'CursorDown 1'
  305.     end
  306. end
  307. 'CursorDown 2'
  308. 'Column -1'
  309. 'Put' format(pr,,4)
  310. 'Column -1'
  311. 'Put "Pearson r:"'
  312. 'CursorDown 1'
  313. 'Put "r sq.:"'
  314. 'Column 1'
  315. 'Put' format(rs,,4)
  316. 'CursorDown 1'
  317. 'Put' format(Sterr,,4)
  318. 'Column -1'
  319. 'Put "Std. Error of r:"'
  320. 'CursorDown 1'
  321. 'Put "N:"'
  322. 'Column 1'
  323. 'Put' NRows
  324. 'CursorDown 1'
  325. 'Put' df
  326. 'Column -1'
  327. 'Put "df:"'
  328. 'CursorDown 1'
  329. 'Put "t:"'
  330. 'Column 1'
  331. 'Put' format(t,,6)
  332. 'CursorDown 1'
  333. 'Put' format(P,,6)
  334. 'Column -1'
  335. 'Put "P(T<=t) one-tail:"'
  336. 'CursorDown 1'
  337. 'Put "T-Critical (95%):"'
  338. 'CursorDown 1'
  339. 'Put "T-Critical (99%):"'
  340. 'CursorDown 1'
  341. 'Put "P(T<=t) two-tail:"'
  342. 'CursorDown 1'
  343. 'Put "T-Critical (95%):"'
  344. 'CursorDown 1'
  345. 'Put "T-Critical (99%):"'
  346. 'CursorUp 4'
  347. 'Column 1'
  348. 'Put' format(TCRIT1,,4)
  349. 'CursorDown 1'
  350. 'Put' format(TCRIT2,,4)
  351. 'CursorDown 1'
  352. 'Put' format(PT,,6)
  353. 'CursorDown 1'
  354. 'Put' format(TCRIT3,,4)
  355. 'CursorDown 1'
  356. 'Put' format(TCRIT4,,4)
  357. 'Column -1'
  358. 'CursorDown 1'
  359. If PCORR~=0 then DO
  360. 'Put "Test Against True Correlation:"'
  361. 'CursorDown 1'
  362. 'Put "True Correlation:"'
  363. 'Column 1'
  364. 'Put' PCORR
  365. 'Column -1'
  366. 'CursorDown 1'
  367. 'Put "z:"'
  368. 'CursorDown 1'
  369. 'Put "Prob (z to Left tail):"'
  370. 'CursorDown 1'
  371. 'Put "Prob (z to Right tail):"'
  372. 'CursorDown 1'
  373. 'Put "Density Function:"'
  374. 'CursorUp 3'
  375. 'Column 1'
  376. 'Put' format(zed,,4)
  377. 'CursorDown 1'
  378. 'Put' format(PZ,,8)
  379. 'CursorDown 1'
  380. 'Put' format(Q,,8)
  381. 'CursorDown 1'
  382. 'Put' format(PDF,,6)
  383. 'Column -1'
  384. 'CursorDown 1'
  385. 'Put "Z Critical One-tail(95%):"'
  386. 'Column 1'
  387. 'Put' 1.65
  388. 'CursorDown 1'
  389. 'Put' 2.33
  390. 'Column -1'
  391. 'Put "Z Critical One-tail(99%):"'
  392. 'CursorDown 1'
  393. 'Put "Z Critical Two-tail(95%):"'
  394. 'Column 1'
  395. 'Put' 1.96
  396. 'CursorDown 1'
  397. 'Put' 2.58
  398. 'Column -1'
  399. 'Put "Z Critical Two-tail(99%):"'
  400. 'CursorDown 1'
  401. End
  402. 'Put "Confidence Intervals"'
  403. 'CursorDown 1'
  404. 'Put "95% (+1.96):"'
  405. 'CursorDown 1'
  406. 'Put "95% (-1.96):"'
  407. 'CursorDown 1'
  408. 'Put "99% (+2.58):"'
  409. 'CursorDown 1'
  410. 'Put "99% (-2.58):"'
  411. 'CursorUp 3'
  412. 'Column 1'
  413. 'Put' format(CIP95,,4)
  414. 'CursorDown 1'
  415. 'Put' format(CIN95,,4)
  416. 'CursorDown 1'
  417. 'Put' format(CIP99,,4)
  418. 'CursorDown 1'
  419. 'Put' format(CIN99,,4)
  420. 'Refresh 1'
  421. 'Refresh 2'
  422. /*'Message' "Finished"*/
  423. /*indicate the main script is finished*/
  424. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  425. result=writeln(6Info, DisplayMsg)
  426. if result~=0 then do
  427.     /*Wait 3 seconds*/
  428.     CALL DELAY(150)
  429.     /* close window*/
  430.     result=close(6Info)
  431. end
  432. 'DEFPUBSCREEN("Workbench")'
  433. exit
  434.  
  435. /* Procedures */
  436.  
  437. cellrow: procedure
  438. do
  439.     parse arg cell
  440.     do charpos=2 to length(cell)
  441.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  442.     end
  443.     return 0
  444. end
  445. Return
  446.  
  447. cellcol: procedure
  448. do
  449.     parse arg cell
  450.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  451.     cell=upper(cell)
  452.     len=length(cell)
  453.     val=0
  454. do charpos=1 to len
  455.     if datatype(substr(cell,charpos,1),n) then
  456.     do cell=reverse(substr(cell,1,charpos-1))
  457.     do x=1 to length(cell)
  458.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  459.     end
  460.     return val
  461.     end
  462.     end
  463.     return 0
  464. end
  465. Return
  466. /* It is important to put the exposed array at the end of the next line */
  467. Sort: procedure expose NCols NRows data.
  468. do x=1 to NCols
  469. L=(xtoy(2,int(log(NRows-1)/log(2))))-1
  470.     Do Until L<1
  471.     L=trunc(int(L/2))
  472.     Do J=1 to L
  473.             Do K=J+L To NRows-1 By L
  474.             I=K
  475.             dumdat=data.x.I
  476.             Do while I>L
  477.                 y=I-L
  478.                 If data.x.y ~> dumdat then Leave
  479.                 data.x.I=data.x.y
  480.                 I=I-L
  481.             End
  482.             data.x.I=dumdat
  483.             End
  484.         End
  485.     End
  486. End
  487. Return
  488.  
  489. syntax:
  490.      if arg(1)='FAIL' then do
  491.         'Message "Library is unavailable."'
  492.         'DEFPUBSCREEN("Workbench")'
  493.         exit
  494.         end
  495.     'DEFPUBSCREEN("Workbench")'
  496.     exit
  497.  
  498. Format:  procedure
  499.  
  500.      arg number, before, after
  501.      CallLine = SIGL
  502.      if ~datatype(CallLine, 'N') then CallLine = '??'
  503.  
  504.      /* Make sure we have a number as first (required) argument    */
  505.      if ~datatype(number, 'N') then do
  506.         if number = '' then
  507.            fc = 17     /* Wrong number of arguments           */
  508.         else
  509.            fc = 47     /* Arithmetic conversion error             */
  510.         signal FormatSyntaxError
  511.      end
  512.      num = number + 0
  513.      if before = '' & after = '' then
  514.         return num
  515.      else do
  516.         parse var num integer '.' fraction
  517.         if before = '' then before = length(integer)
  518.         if after = '' then after = length(fraction)
  519.         if ~datatype(before, N) | ~datatype(after, N) then
  520.            do fc = 18
  521.            signal FormatSyntaxError
  522.        end
  523.         if before < length(integer) then do
  524.            fc = 18
  525.            signal FormatSyntaxError
  526.         end
  527.         if after ~= length(fraction) then do
  528.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  529.         if integer<1&integer>-1 then integer=integer
  530.            else integer = integer + (fraction % 1)
  531.            fraction = substr(fraction, 3)
  532.         end
  533.         if fraction >= 0 then
  534.            return right(integer, before)'.'fraction
  535.         else
  536.            return right(integer, before)
  537.      end
  538.  
  539.  FormatSyntaxError:
  540.         if show('F', STDERR) then
  541.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  542.         else
  543.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  544.         'Message' mess
  545.         parse source Func .
  546.         if Func = 'FUNCTION' then do
  547.         'DEFPUBSCREEN("Workbench")'
  548.            exit "Err"
  549.         end
  550.         else do
  551.         'DEFPUBSCREEN("Workbench")'
  552.            exit 10
  553.         end
  554.  
  555. PROBT: PROCEDURE
  556.  
  557.     PARSE ARG TA,K1
  558.     A=.36338023
  559.     W=ATAN(TA/SQRT(K1))
  560.     S=SIN(W)
  561.     C=COS(W)
  562.     L=K1-2*INT(K1/2)
  563.     IF L=0 THEN DO
  564.         T1=S
  565.         IF K1=2 THEN DO
  566.         PO=.5*(1+T1)
  567.         RETURN PO
  568.         END
  569.         T2=S
  570.         J1=-1
  571.         J2=0
  572.         K2=(K1-2)/2
  573.     END
  574.     ELSE DO
  575.         T1=W
  576.         IF K1=1 THEN DO
  577.             T1=T1*(1-A*L)
  578.             PO=.5*(1+T1)
  579.             RETURN PO
  580.         END
  581.         T2=S*C
  582.         T1=T1+T2
  583.         IF K1=3 THEN DO
  584.             T1=T1*(1-A*L)
  585.             PO=.5*(1+T1)
  586.             RETURN PO
  587.         END
  588.         J1=0
  589.         J2=1
  590.         K2=(K1-3)/2
  591.     END
  592.     DO I=1 TO K2
  593.         J1=J1+2
  594.         J2=J2+2
  595.         T2=T2*C*C*J1/J2
  596.         T1=T1+T2
  597.     END
  598.     T1=T1*(1-A*L)
  599.     PO=.5*(1+T1)
  600. RETURN PO
  601.  
  602. TCRITICAL: PROCEDURE
  603.  
  604.     PARSE ARG PO,K1
  605.     AO=.0001
  606.     E=.005
  607.     E2=E+E
  608.     A=2*PO-1
  609.     IF K1=1 THEN DO
  610.         TO=TAN(1.57079633*A)
  611.         RETURN TO
  612.     END
  613.     IF K1=2 THEN DO
  614.         SN=SIGN(2/(1-A*A))
  615.         IF SN=-1 THEN TO=-1*(A*SQRT(ABS(2/(1-A*A))))
  616.         ELSE TO=A*SQRT(2/(1-A*A))
  617.         RETURN TO
  618.     END
  619.     P1=PO
  620.     Z=NORMALPP(PO)
  621.     W=Z*(1+2/(1+8*K1))
  622.     T3=K1*(EXP(W*W/K1)-1)
  623.     SELECT
  624.         WHEN Z<0 THEN TT=-1
  625.         WHEN Z=0 THEN TT=0
  626.     OTHERWISE TT=1
  627.     END
  628.     SN=SIGN(T3)
  629.     IF SN=-1 THEN T3=-1*(TT*SQRT(ABS(T3)))
  630.     ELSE T3=TT*SQRT(T3)
  631.     /* LABELA:
  632.     TO=T3
  633.     PO=PROBT(TO,K1)
  634.     F1=PO-P1
  635.     TO=T3+E
  636.     PO=PROBT(TO,K1)
  637.     F2=PO
  638.     TO=T3-E
  639.     PO=PROBT(TO,K1)
  640.     F2=F2-PO
  641.     F2=F2/E2
  642.     T4=T3-F1/F2
  643.     IF ABS(T4-T3)>AO THEN DO
  644.         T3=T4
  645.         SIGNAL 'LABELA'
  646.     END
  647.     */
  648.     T4=0
  649.     DO FOREVER
  650.         TO=T3
  651.         PO=PROBT(TO,K1)
  652.         F1=PO-P1
  653.         TO=T3+E
  654.         PO=PROBT(TO,K1)
  655.         F2=PO
  656.         TO=T3-E
  657.         PO=PROBT(TO,K1)
  658.         F2=F2-PO
  659.         F2=F2/E2
  660.         T4=T3-F1/F2
  661.         IF ABS(T4-T3)>AO THEN T3=T4
  662.         ELSE LEAVE
  663.     END
  664.     TO=T4
  665.     
  666. RETURN TO
  667.  
  668. LOGGAMMA: PROCEDURE
  669.  
  670.     ARG XA
  671.     C1=76.18009173
  672.     C2=-86.50532033
  673.     C3=24.01409822
  674.     C4=-1.231739516
  675.     C5=.001208580
  676.     C6=-.000005364
  677.     C7=2.506628275
  678.     X1=XA-1
  679.     W=X1+5.5
  680.     W=(X1+.5)*LN(W)-W
  681.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  682.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  683.     L=W+LN(C7*S)
  684. RETURN L
  685.  
  686. NORMALPP: PROCEDURE
  687.  
  688.     ARG P0
  689.     A1=2.515517
  690.     A2=.802853
  691.     A3=.010328
  692.     A4=1.432788
  693.     A5=.189269
  694.     A6=.001308
  695.     Q=.5-ABS(P0-.5)
  696.     SN=SIGN(-2*LN(Q))
  697.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  698.     ELSE W=SQRT(-2*LN(Q))
  699.     W1=A1+W*(A2+A3*W)
  700.     W2=1+W*(A4+W*(A5+A6*W))
  701.     ZZ=W-W1/W2
  702.     SELECT
  703.         WHEN (P0-.5)<0 THEN TT=-1
  704.         WHEN (P0-.5)=0 THEN TT=0
  705.         otherwise TT=1
  706.     END
  707.     ZZ=ZZ*TT
  708. RETURN ZZ
  709.  
  710. NPROB: Procedure Expose PZ Q PDF
  711.  
  712.     arg Z
  713.     A0=0.5
  714.     A1=0.398942280444
  715.     A2=0.399903438504
  716.     A3=5.75885480458
  717.     A4=29.8213557808
  718.     A5=2.62433121679
  719.     A6=48.6959930692
  720.     A7=5.92885724438
  721.     B0=0.398942280385
  722.     B1=3.8052E-8
  723.     B2=1.00000615302
  724.     B3=3.98064794E-4
  725.     B4=1.98615381364
  726.     B5=0.151679116635
  727.     B6=5.29330324926
  728.     B7=4.8385912808
  729.     B8=15.1508972451
  730.     B9=0.742380924027
  731.     B10=30.789933034
  732.     B11=3.99019417011
  733.     ZABS = ABS(Z)
  734.     Y = A0*Z*Z
  735.     PDF = EXP(-Y)*B0
  736. SELECT
  737.     WHEN (Z>=-1.28) & (Z<=1.28) then DO /*Z BETWEEN -1.28 AND +1.28*/
  738.         Q = A0-ZABS*(A1-A2*Y/(Y+A3-A4/(Y+A5+A6/(Y+A7))))
  739.         IF (Z<0) then do
  740.             PZ = Q
  741.             Q = 1.0-PZ
  742.         End
  743.         Else PZ = 1.0-Q
  744.         RETURN
  745.     End
  746.     WHEN (Z>1.28) & (Z<=12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
  747.         Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
  748.         IF(Z<0) then Do
  749.             PZ = Q
  750.             Q = 1.0-PZ
  751.         End
  752.         Else PZ = 1.0-Q
  753.         RETURN
  754.     End
  755.     WHEN (Z<-1.28) & (Z>=-12.7) then do /*ZABS BETWEEN 1.28 AND 12.7 */
  756.         Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/(ZABS+B9+B10/(ZABS+B11))))))
  757.         IF(Z<0) then Do
  758.             PZ = Q
  759.             Q = 1.0-PZ
  760.         End
  761.         Else PZ = 1.0-Q
  762.         RETURN
  763.     End
  764.     OTHERWISE /*Z FAR OUT IN TAIL*/
  765.         Q = 0
  766.         PDF = 0
  767.         IF(Z<0) then do
  768.             PZ = Q
  769.             Q = 1.0-PZ
  770.         End
  771.         Else PZ = 1.0
  772.         RETURN
  773. End
  774. Return